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ABSTRACT 



Context. The gravitational wave (GW) background in the range 0.01 - 30 mHz has been assumed to be dominated by unresolved 

radiation from double white dwarf binaries (DWDs). Recent investigations indicate that, at short periods, a number of DWDs should 

be resolvable sources of GW. 

Aims. To characterize the GW signal which would be detected by LISA from DWDs in the Galaxy. 

Methods. We have constructed a Galactic model in which we consider distinct contributions from the bulge, thin disc, thick disc, and 

halo, and subsequently executed a population synthesis approach to determine the birth rates, numbers, and period distributions of 

DWDs within each component. 

Results. In the Galaxy as a whole, our model predicts the current birth rate of DWDs to be 3.21 x 10"^ yr ', the local density to be 

2.2x lO^^'pc"^ and the total number to be 2.76 x 10*. Assuming SNIa are formed from the merger of two CO white dwarfs, the SNIa 

rate should be 0.0013 yr"'. The frequency spectra of DWD strain amplitude and number distribution are presented as a function of 

galactic component, DWD type, formation channel, and metallicity 

Conclusions. We confirm that CO+He and He+He white dwarf (WD) pairs should dominate the GW signal at very high frequencies 

(log/Hz"' > -2.3), while CO+CO and ONeMg WD pairs have a dominant contribution at log/Hz"' < -2.3. Formation channels 

involving two common-envelope (CE) phases or a stable Roche lobe overflow phase followed by a CE phase dominate the production 

of DWDs detectable by LISA at log /Hz"' > -4.5. DWDs with the shortest orbital periods will come from the CE+CE channel. The 

Exposed Core plus CE channel is a minor channel. A number of resolved DWDs would be detected, making up 0.012% of the total 

number of DWDs in the Galaxy. The majority of these would be CO+He and He+He pairs formed through the CE+CE channel. 
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1. Introduction 

Substantial efforts have been made to measure gravitational 
waves that, according to Einstein, should permeate the whole 
Universe (Press & Thorne 1972; Tvs on & Giffa rd 1978). Today, 
■ four major ground-based detectors are in operation (LIGO, 
VIRGO, GEO600, and TAMA300), but no direct evidence for 
the existence of gravitational waves (GW) has been established. 

Following the proposed launch of LISA (Laser 
' Interferometer Space Antennae), it is anticipated that the 
most easily detected waves will originate from close compact 
binary stars. In addition to providing a direct test of general rela- 
tivity, the observations of GW from such systems will contribute 
to understanding close binary star evolution, the distribution of 
X-ray sources, y-ray bursts and supernovae explosions, and may 
also improve knowledge of galactic structure. 

As a result of recent and cuiTent surveys, increasing num- 
bers of candidate compact binaries are being found, including 
double neutron stars (NS) and double white-dwarf (DWD) sys- 
tems. AMCVn variables are thought to be semi-detached ultra- 
compact binary stars systems in which a white dwarf accretes 
from a degenerate (or semi-degenerate) companion with an or- 
bital period of less than 80 minutes. CuiTently the nu mber of 
knowi i AMCVn sy stems is approximately 20 (Ram sav et alJ 
I2007t iRoelofs et al]|2007l) . while those of detached DWD and 
NS systems with known orbital periods are 24 and 8 respectively. 



Theoretical studies of the dominant GW s ignal in the 
Galaxy have been pursued for several decades. ' Mironovskiil 
(1966) calculated the spectral density and total flux due to 
GW from W UMa binaries and predicted a very low flux with 
a spectrum peaked at a period of about 0.17d (0.07 mHz) 
(A- and W-type WUMa binaries have orbital periods in the 
rang e 0.4 - 0.8 d and 0.22 - 0.4 d, respectively). Subsequently 
Eva ns et al] ([1987 ) demonstrated that degenerate dwarf binaries 
should be detectable sources of gravitational waves. More re- 
cent work indicates that DWD binaries a re cruc i al GW sources 

19901 iHils & Bendeil 



in th e 0.1-1 00 mHz range dHiIsetaL. 

2000'; 'Nelemans et alj |2001bl |2004'; 'Farmer & Phinney' |200l 



Willems et al. 200'O'he results of Ruiter et al. (2009) show the 
halo signal would have a negligible effect on the detection of 
LISA sources in comparison with that of the disc and bulge. 

Such population synthesis models are sensitive to the 
adopted inputs, including the fast stellar evolution models, the 
Galactic model, including mass and metallicity distribution, ini- 
tial mass function, and star formation rate, and various parame- 
ters governing binary star evolution, such as mass loss and ex- 
change in binary stars. It is important to ensure that previous 
results are robust by reproducing them in independent calcu- 
lations, and to identify the major areas of uncertainty in these 
models. It is also important to validate a new code by testing it 
against established benchmarks, before using it to explore new 
physics. This paper addresses both of these functions. 
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To discuss the GW background due to a population of galac- 
tic DWDs, we have executed a simulation based on a popula- 
tion synthesis approach for DD binaries which uses a multi- 
component model of the Galaxy. We do not include the contri- 
bution of DD's containing an NS or black hole (BH) since the 
formation and evolution of these await further theoretical and 
observational constraint. Other parameters of our stellar evolu- 
tion model are taken from the best model of Han ( 1998). 

In addition to the GW signal, our model provides an esti- 
mate of the current birth rate, total numbers and properties (e.g. 
the primary mass and orbital period) of DWDs which may be 
compared with observation. 

In §2 we describe the model for generating a GW signal from 
a binary, the stellar evolution model for generating a population 
of compact binaries, and the model for generating their galactic 
distribution. In §3, we show the main results from these mod- 
els in terms of the overall Galactic GW signal, and in terms of 
the contribution from various components of the DWD popula- 
tion. Section 4 compares these results with previous work, and 
conclusions are drawn in §5. 



(iDouglass & BraginskvllT979t iNelemans et"ani2001bl) . So from 
EqslT]|6] and Kepler's law a^/P^^^ - G(mi + m2)/4-n^, we can 
write li(n, e) in the n* harmonic as 



h(n,e)^ 1.0x10" 



X 



g(n,eW"lM''^' 



M\' (Porb\ 

Mel \ h / 



-2/3 /^\-l (7) 

kpc/ 



Equation |7] gives the strain amplitude in the «* harmonic with 
a frequency /^^ - n/Porb, and where M = /r'^^M^^^ is the so- 
called chirp mass. Especially for a circular orbit, i.e. e — Q, n - 
2, we have 



/i(n,e) = /!(2,0) = 5.0x 10 

\5/3 ,D 2/3 



22 



MqI I h / \kpc/ 



(8) 



2.2. Stellar evolution models 



2. Simulations 

2.1. The GW signal from a binary star 

The power of the GW radiated from two point masses in the n* 
harmonic is given by 

32/GV'Mn 
LGw(n,e)=— — — —\g(n,e) (1) 

where G is Newton's gravitational constant, c is the speed of 
light, a is their orbital separation, e is the eccentricity of the or- 
bit, M and fi are the total and reduced masses of the system, and 
g{n, e) is a coefficient expressing the relative pow er in each har- 
monic and is given bv iPeters & Mathewsl (Il963h . If mi and m2 
are the masses of the components, then 



M - m\ + m2, i-i — m\m2lM. 
The GW energy flux density incident on the Earth is 



icw _ c^ 
4^ ~ 16;rG 



< h\+h\ >, 



(2) 



(3) 



where d is the distance, h+ and hx are dimensionless functions 
for the amplitude of the wave i n two orthogonal polarizations, 
and dots denote time der ivatives (llsaacsonll968l : lPress & Thornel 
119721: Evans et al.lll987l) . Assuming a binary system in a circular 
orbit of period f orb at a distance d, h+ and h^ are given by 



h. 



c4 d 



2(1 +cos^i)i7:fcwMf^^ficos(27TfGwt), (4) 



(55/3 J 
hx = +^^-74cos/(7r/GwM)2/Vsin(2;r/GwO, 



(5) 



where /gw = 2/f orb is the frequency of the emitted GW, and / is 
the inclin ation angle (; = 90° indic ates a binary which is viewed 
edge-on) (iLandau & Lifshitzlll975 !). 

The so-called strain amplitude h, defined as h^ - j [h\ ^^y. -H 
hxmaxl' describes the strength of the signal to be measured by 
GW detectors. From Eqs.[3]|4]and|5] we obtain 



4GLgw 



1/2 



(6) 



A binary star population synthesis calculation requires the pre- 
dicted evolution of many thousands of binary stars, only a frac- 
tion of which yield systems of interest 02.2.6I I. We have ob- 



tained this predicted evolution using the fas t binary star evo- 

aD d2000H2002h and 



lution (BSE) code descri bed bv iHurlev et aL 

previously referred to bv iHanl (Il998l) . However, we have modi- 
fied the code in a number of ways, including several of the ap- 
proximation formulae and the treatment of the common enve- 
lope ejection mechanism. This section describes these changes, 
including the approximations used to simulate changes in mass 
and angular momentum due to the interaction of binary star com- 
ponents. 

2.2.1. Common-envelope (CE) evolution 

RLOF occurs in binary systems when one star fills its Roche 
lobe either by evolutionary expansion or by orbital shrinkage. 
The Roche radius of the primary is given by 



a 



0.49q 



2/3 
1,2 



.2/3 



l/3^ 



0.6^t,2+ln(l+?i2) 



(9) 



where a is more generally the s emi-major axis of the orbit and 
'71,2 = 'Wi,2/'W2,i (Eggleton 1983). Studiesof the critical value for 
stable RLOF have shown that q^ is a function of primary mass 
nil, its core mass mic, and mass-transfer efficiency of the donor. 
We adopt 



qc 



1.67 -xH- 2 



mi 



/2.13 



(10) 



where x - 0.3 is the exponent o f the mass-ra(3ius relation at con- 
stant luminosity for giant stars (" Hurley et al.ll200(il2002h . 

The calculation of the properties (e.g. orbital separation) of 
a binary after CE ejection in our model is based on the varia- 
tion of angular momentum (y-mechanism). In the case of non- 
conservative mass transfer, the angular momentum loss of a bi- 
nary system would be described by the decrease of primary mass 
times a y factor (£aczvnski & Zidlkowski 1967j Nelemans et aP 
I2OOOI) : 

Ji-Jf mi-mic ^,,^ 

7 ' (11) 



Ji 



M 



where Ji is the orbital angular momentum of the pre-mass trans- 
fer binary, and 7f is the final orbital angular momentum after CE 
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ejection. If initial and final orbits are assumed circular, the frac- 
tion of angular momentum lost during the mass transfer, 7i - 7f , 
becomes 



Ji - Ji 



niynii 2 

a- O): - 

M ' 



micm2 
(mic + OT2) 



a^cDf 



(12) 



where w, and ojf are the circular frequency of the binary before 
and after CE ejection. Combining fTTI and fT2l with Kepler's law, 
we have the ratio of final to initial orbital separation 



Of _/»Jl_\ / OTlc +'W2 \/ 



mi — mic 



(13) 



M n M / 

iNelemans & Tout! (l2005l) investigated the mass-transfer phase of 
the progenitors of white dwarfs in binaries employing the 7- 
mechanism based on 10 observed systems and suggested the 
value of 7 in the range of 1.4 ~ 1.7. We here adopt y - 1.5, 
following the recommendation by INeleman s et al. (2000). 

In order to eject a binary star common envelope, conserva- 
tion of energy requires that the envelope binding energy, includ- 
ing gravitational binding a nd recombination ener gies, must be 
equal to the orbital energy (Webbink"l984", '2008"). The ratio of 
final to initial orbital separation can be expressed as 



Of 

a; 



mi 



1 +2 



1 



1 



XefiaiX mi -OTic 



/icE^nRhi /ip^Li Gm 



m2 



(14) 
where the coefficients are given bv lWebbinkI (l2008h . 

Both CE ejection mechanisms can reproduce observation s 
(INelemans et al.ll2000t INelemans & Tbuil200l IWebbinkll2008l) . 
We here adopt the y-algorithm in order to reduce the number of 
free parameters in the model. 

We note that there is a major difference between the a- and 
■y- algorithms. In both, the first CE phase can lead to the pro- 
duction of a low-mass binary with q x I. However, the a- al- 
gorithm requires a significant spiral-in stage in order to eject the 
envelope whilst the y- mechanism does not, implying that the 
orbit al separation can be larger after CE ejection in the latter 
case ( INelemans et aliUo OO; Webbin lc..2008.) . 



2.2.2. Stellar wind and accretion 

We assume that the existence of a close companion will in- 
crease the rate of mass-loss (m) due to a stellar wind from a 
star in a binary system. The tidal enhancement of ihi is mod- 
elled by Reimersl ( 1975h formula with an extra tidal term by 
lTout&Eggletonl(il988l) : 



nil 



-4x 10 



mi 



1 + Bx |min(— , -) 

ru 2 



6\ 



yr 



(15) 



where ri, li and mi are stellar radius, lum inosity and mass in 
solar units, respectively. We take B - 1000 (lHanlll998h . 

Some of the mass lost in the wind may be accreted 
by the c ompanion. The classi cal accretion rate formula is 
given by 'Bondi & Hoyy (1 19441) . and the mean accretion rate 
(.Boflin & Jorisse n 198Cis 



-1 



m2 



vr^ 



1 



Gm2 \ j6w 

vl I 2fl2 (1 + y2)3/2 



mi. 



(16) 



where 1 < By, < 2 is an arbitrary parameter (J3s^ 
(iBoffin et al.l Il993l) ). Vw is the wind velocity, v^ - 



orb 



3/2, 



and Vorb = (GM/a) ' is the orbital velocity. Here, we adopt 



Vw = 20kms '. The angular momentum of a star will change if 
the star loses or gains mass. We assume that the specific angular 
momentum in the wind is proportional to the angular velocity 
at the surface of the mass-losing star, and that th e wind angu- 
lar momentum is transferred with 100% efficiency dHurlev et al.l 
2002). 



2.2.3. Gravitational radiation and magnetic braking 

A close compact binary system driven by gravitational radiation 
would eventually undergo a mass transfer phase, ultimately lead- 
ing to coalescence. The orbital changes due to gra vitational ra- 
diatio n from two point masses are predicted to be dHurlev et al.l 
I2OOI : 



a 
2a 



•'gr 

■'orb 



8.315x10" 



ioim_ 

Mo 



m2 

Mo 



M 

Mo 



1-^F ] 

(l_e2)5/2^ 



(17) 



yr 



8.315x10 



-10/^ 

Mr 



m2 

Mo 



M 

Mo 



29 , 121 2 \ 

6 "^ 96 ^ 

(l_e2)5/2^ 



(18) 



yr 



Gravitational radiation could explain the formation of cata- 
clysmic variables (CVs) with orbital periods less than 3h, while 
magnetic braking of the tidally coupled primary by its own mag- 
netic wind would account for orbit al angular-momentum loss 
from CVs with periods up to 10 h riFaulknedl971 ^JZangrilli et al.l 
1997). We use the formula for the rate of angu lar-momentum 
loss due to magnetic braking derived bv lRappaport et al.l (119831) 
and .Skumanich (.1972.) : 



■,me 
m 



j^,^-5.S3xlO-^''-^^\—^\ MoR^yr 

Koyr 



(19) 



where r, menv and m are the radius, envelope mass and mass of 
a star with a convective envelope, and cospin is the spin angular 
velocity of the star. 

2.2.4. Tidal interaction 

Observations indicate that stellar rotation in a clo se binary sys- 
tem tends to synchronize with the orbital motion (lLevatolll974t 
IStrassmeienfl996) . The two tidal dissipation processes, i.e. tur- 
bulent dissipation of the equilibrium tide and radiative damping 
of the dynamical tide, are able to account for the observed orbital 
circularization of close binary and components rotating in syn- 
chronism with orbital motion; see Zahn (2 005b for a review. The 
calcul ation of tidal interaction in our model follows JHurlev et aP 
(I2OO2I) and adopts an equilibrium t ide for the stars with a con- 
vective envelope ( Rasio et al . '1996*) and a dynamical tide for the 
stars with a radiative envelope (Zahn 1975, 1977). 

In addition, we note that a hydrodynamical mechanism may 
account for the orbital circularization as it is concomitant with 
the hydrodynamical spin-down of the com ponents both in early- 
type and late-type detached close binaries ( Tassoullll988l) . Since 
this mechanism would be more effective than the tidal interac- 
tion at circularizing the orbit of a binary, we adopt the tidal dis- 
sipation as an upper limit to the circularization timescale. 



S. Yu & C. S. JefFery: Gravitational Waves from Double White Dwarfs 



The tidal synchronisation time is then 






5/7 



m 



yi- 



(20) 



for a tide raised on a white dwarf secondary of mass mo 
I II 1 

(ICampbeUII 19841) . where r2, h and q2 represent its radius, lumi- 
nosity and mass ratio respectively, fsynch would be a lower limit if 
tidally-induced non-radial oscillations were included. The orbit 
of DWDs should be circularized akeady, and a circularization 
time-scale is not relevant. However, the synchronization time- 
scale is essential as the companion may be spun up by mass 
transfer dHurlev et al..i2002) . 

2.2.5. Mass transfer in compact binaries 

As a result of gravitational radiation, mass transfer or direct 
merger between two compact binaries is inevitable if the com- 
pact object with lower mass fills its Roche lobe. The ana- 
lytic mass -radius relation for zero-temperature white dwarfs, is 
(lNauenberg.1 972) 



r = O.O1125/?0(ot' 



-2/3 



,/2/3y/2 



(21) 



where m' - m/1.433M0 is the mass in terms of the 
Chandrasekhar mass. Combining with Eq. |9l we deduce an equi- 
librium mass transfer rate 



'«2 



C2^Edd/0i-l - 2m2/TGR 
^ad - Ci - C2Cf>L2/(f>il 



(22) 



where Lem = Anr^cglK is the Eddington luminosity, k - 
0.2(1 + X)cm^g"' is the opacity of the accreted gas, assumed 
here to be due to electron scattering, X is the hydrogen mass 
fraction, g is the g ravitational acceleratio n, and other coefficients 
are as given by iHan & WebbinkI (Il999l) . tgr is the time scale 
for orbital angular-mome ntum loss due to gravitational radiation 
dLandau & Lifshitzll 19581) with 



Tgr = 



5 c' 



32G^ mim2M' 



(23) 



and Ci, C2 ^ad and 0li as given by Han & Webbink ( 1999). 

For stable mass transfer, .Han & Webbink (.1999.) give the 
critical mass ratio 



q = — < q, ^ 0.7 - Q.l(mi/Me). 
m\ 



(24) 



Here, m2 is donor. lfq> q^, the mass transfer becomes dynami- 
cal ly unstable, even t ually c ausing a runaway merger. 

iNelemans et al.l (12001 ah give an alternative value q^ ~ 5/6 -H 
^(m2)/2, ^ = dlnr/dlnm, using the size of the R oche lobe and 
the mass transfer rate de rived bv lPaczviiskil (Il967l) . We adopt the 
iHan & WebbinkI (Il999h approximation. 

Finally, we give a formula for the Eddin gton accretion rate 
(ICameron & Mock|[T967HHurlev et al.ll200l 



MEdd 



2.08 X 10^^(1 -HX)"ViM0yr"', 



(25) 



which places a significant limit on the amount of mass accreted 
by compact objects. 

We set this limit only for the accretion of compact ob- 
jects. However, whether the Eddington limit should be applied 
or not needs to be further investigated, since the luminosity 
generated by accretion in excess of this limit might be carried 
away from the system via a strong wind, or transfered to a 



disc. Super-Eddington accretion rates may be important for the 
formation of low-mass X-ray binaries and millisecond pulsars 
(Webbink & Kalogera 1997) and also for modeling X-ray emis- 
sion from quasars due to accretio n from a disc onto a fast rotating 
black hole (.Beloborodovll 19981) . I mposing the limit diminishes 
the birth rate of type la supernovae (lLivioll200u) . 



2.2.6. Tine formation clianneis for compact binaries 

One can find a detailed discussion of the formation channels of 

compact binaries in Han (1998). We here only summarize the 

three main evolution channels which produce the majority of 

compact binaries : 

I: Stable RLOF + CE ejection; 

II: CE ejection + CE ejection; 

III: Exposed core + CE ejection. 

In channel I, a binary with mass ratio q - mi 1 1112 less than 
some critical value q^ will experience dynamically stable mass 
transfer if the primary fills its Roche lobe while the star is in the 
Hertzsprung gap or on the red giant branch. The primary will 
become a compact object and the orbital separation will change 
as 

- d In a = 2d In m2 -H 2aRL0Fd In mi + d \vi(m\ + OT2) (26) 



for stable Roche- 
Here, we take 



where ffRLOF is the mass-transfer effi ciency 

lobe overflow (RLOF) (Hanetal. '199 51) 

ffRLOF = 0.5 (Paczvhski &Zi61kowski 119671 : iRefsdal et al.l 

1974). Subsequently, if the secondary fills its Roche lobe while it 

is in the Hertzsprung gap or on the red giant branch, then RLOF 

will occur 

If the adiabatic response of the radius of the mass donor is 
less than the change of its Roche lobe radius with respect to a 
change of mass, i.e. ( ^,'"5'°"'" ) < ( 1'"^^-°" ) , mass trans- 
fer will be unstable and a common envelope (CE) will form. 
Interaction (friction) between the compact cores and the CE will 
convert orbital energy into kinetic energy, heating and expand- 
ing the CE. If the energy conversion mechanism is sufficiently 
efficient, the CE will be expelled and a compact binary with a 
short orbital period will result. 

Channel II demands that the binary has q > q^ so that, if 
the primary fills its Roche lobe while in the Hertzsprung gap or 
on the red giant branch, unstable mass transfer will lead to the 
formation of a CE. Similarly, if sufficient orbital energy can be 
extracted, the CE will be ejected to leave a binary containing 
a compact object. Following evolution of the secondary away 
from the main sequence, a second CE could form and be ejected 
to leave a compact binary. 

Channel III represents a variation of channels I and II in 
which the envelope of a massive primary is removed by a stellar 
wind rather than a first CE ejection. CE ejection following evo- 
lution of the secondary may also give rise to a compact binary. 

I n sum mary, CE plus CE represents the major channel. Up to 
74% (lHan.1998,) . or 67% in our model (see Table©, of compact 
binaries may be generated through this channel depending on the 
choice of parameters. This may also be the most efficient channel 
to produce compact binaries with very short orbital periods. The 
stable RLOF plus CE ch annel can account for up to 50% of com- 
pact binaries (lHanlll998,) . or 30% in our model. Exposed core 
plus CE is a minor channel which could make a significant con- 
tribution when a binary evolves with a strong tidally-enhanced 
stellar wind. Both the results of Han ( 1998) and our own show 
this channel to be negligible for compact binaries with orbital 
periods shorter than 70 hr 
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Table 1. Density laws and associated parameters, r is the spher- 
ical radius from the center of the Galaxy and ro is bulge scale 
length; R and z are the natural cylindrical coordinates of the ax- 
isymmetric disc, h^ is the scale length of the disc, /i, is the scale 
height of the thin disc, hi is the scale height of the thick disc; a is 
the radius of the halo and aq is a constant; pc is the central mass 
density. 



density law 



constants 
(kpc) 



(Mopc-^) 



Bulge 



-(r/ro^ 



ro = 0.5 f\ = 12.73 



Thin disc e-'^"'''sech^(-z/hj hf, = 2.5 

h, = 0.352 



Thick disc 



Halo 



-Rlhn„-JH 



[(l + (t)')]' 



JiR = 2.5 



1.881 



0.0286 



K 



aa 



1.158 

= 2.7 



0.108 



2.3. Galactic structure 

A detailed model of the Galaxy is important in order to describe 
the overall distribution of white dwarf binary systems, including 
the distance of each from the Sun. It is believed that the cur- 
rent Galaxy mainly comprises the bulge, the thin disc, thick disc 
and the halo. We summarize our approximation for the Galactic 
density distribution in Table [1] 

Figure [T] shows the distributions of distances to the Earth of 
compact binaries in the Galaxy and the bulge plus the thin disc. 
We assume that the position of the sun is /?sun = 8.5 kpc and Zsun 
= 16.5 pc (Eeudenreich 1998). 

(1) We adopt a normal density di stribution for the sphe rical 
bulge with a cut-off radius of 3.5 kpc (JNelemans et al.ll2004l) . 



p,(r) = i^,-(^/^o)= 



4nrl 



Mopc" 



(27) 



where r is the radius from the center of the Galaxy, ro - 0.5 kpc 
is the bulge sca le length, and My, i s the mass of bulge (see ^2.41 
find the value). iRobin et all (12003 ^ suggest that the structure of 
the inner bulge (< 1° from the Galactic center) is not yet well 
constrained observationally. Consequently we here focus on the 
outer bulge and make no allowance for any additional contribu- 
tion to the compact-binary population from the central region. 

(2) A more complicated function is involved in the spatial 
density distribution of stars in the disc. IS ackettI (|l997h p roposed 
a double exponential distribution. Phleps et al] ( l200(]h derived 
three functions for the star density distribution in their model of 
a thin disc plus thick disc (exponential + exponential, hyperbolic 
secant + exponential, and squared hyperbolic secant + exponen- 
tial, respectively) from fits to deep star counts c arried out in the 
Calar Alto Deep Imaging Survey. iRobin et aP (l2003h created a 
complicated function to construct the structure of the Galactic 
disc which is in agreement with Hipparcos results and the ob- 
served rotation curve. 

We here model the thin and thick disc components of the 
Galaxy using a squared hyperbolic secant plus exponential dis- 
tribution expressed as: 



Pd{R,z) 



Ma 
Anh\h 



e-'^/i"<p(z) Mopc 



-3 



(28) 




-5 5 
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Fig. 1. The top panel shows the distribution of stars in the 
R(Galactic plane)-z(height) diagram for our Galactic model at 
age 10 Gyr, using a random sample of ~ 70, 000 stars. The lower 
panel gives the number of DWDs in the Galaxy as a function of 
galactic radius, using a bin size of 0. 1 kpc. The inset panel shows 
the thick disk and halo distribution. We assume there is no immi- 
gration between each Galactic component, i.e. no mass transfer, 
no angular momentum transfer, no collision. 

where R and z are the natural cylindrical coordinates of the ax- 
isymmetric disc, and /i« = 2.5 kpc is the scale length of the disc, 
h = h, for the thin disc, h = h', for the thick disc, and Mj - Mm 
is the mass of the thin disc; M^ - M± is the mass of the thick 
disc (see ^2.4l find the values). p(z) is the distribution in z, with: 



and 



p(z) - sech (-z/h,) (thin disc) 



p(z) = e '^1'^'- (thick disc). 



(29) 



(30) 



where h, - 0.352 kpc is the scale height of the thin disc and 
K - 1.158 kpc is the scale height of the thick disc. We neglect 
the age and mass dependence of the scale height. 

(3) For the halo, we employ a relatively simple density distri- 
bution which is co nsistent with Caldw ell & Ostr iker (.1981) and 
IRobin etalJ(l2003l) :n 



(31) 



Ph(a) = Pc, X 1 + — 

\flO 



where a is the radius of the halo, pch = 0.108 Mopc -^ and 
flo = 2.7 kpc. 



' Alt hough a distribution followin g e"" or a''"* might be a better 
choice dde Vaucouleursi [19531 , Il958h . we have chosen the form given 
here for simplicity. 
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2.4. Population synthesis parameters In the Monte Carlo 
approach 

In order to obtain a sample of compact binaries in the Galaxy, 
we have performed a Monte Carlo simulation in which we need 
five physical inputs: 

(i) We assume that the star formation rate (SFR) in the bulge 
and thin disc is the combination of a main star forming process 
(the first item of the following function) and a minor star forma- 
tion (the second item of the function), 



SFR(0 = lie 



-('-f0)/T 



+ O.12(f-fo)M0yr-',f >fo (32) 



where t is time since the halo was formed. Assuming the current 
age of the Galaxy is 14 Gyr, to - 4 Gyr defines the age of the 
bulge and thin disc to be 10 Gyr and t = 9 Gyr yields a current 
SFR = 4.82 Moyr-i, 1.45 Moyr"' in the bulge and 3.37 Moyr^i 
in the thin disc, on average. The se values are consistent wit h 
ISmith et al.l(ll978h . lTimmes et al.\(T991) and Piehl et al.l(l2006h . 
We assume that SFR(f) = in the bulge and thin disc when 
< f < fo- 

In terms of these assumptions, we infer that the combined 
mass of thin disc and bulge approaches 7.2 x lO"* Mq, shghtly 
higher than the mass of 7.0 x 10'" Mq reported bv .Klvpin et al.l 
(120021) . Within this mass, the bulge contains Mb = 2.0 X 10'° Mq 
and the remaining Mm = 5.2 x 10'" M© is in the thin disc. 

We suppose a burst of star formation, effectively a 6 function, 
happ ened att = Q Gyr f or the halo, and at f = 3 Gyr for the thick 
disc (iRobin et al.ll2003l) and no star formation thereafter. We as- 
sume that the thick disc and the halo attain baryonic masses of 
Mtk = 2.6 X 10^ Mq (5% of thin disc) and Mh = 1.0 x lO"* Mq 
respectively. These numbers are adopted for simplicity solely to 
estimate their contribution to the GW signal. 

(ii) The initial mass function (IMF) can be constrained by the 
local luminosity function, stellar density and potential. We here 
adopted the IMF for the Galacti c compon ents based on the re- 
sults of Robin et al. (2003) and Kroupa et a l. ( 1993) constrained 
by the observations of IWielenetal.1 (11983). Pot^cer (1980) 
and th e Hipparcos mission dCreze et al.lll998i ; iJahreiB & Wieleni 
Il997h . ^ 

Fo r the bulge, we suppose an IMF following iRobin et aP 
(l200l . 



^(in) oc m 



-2.35 



m > 0.7Mp 



(33) 



where m is the primary mass and ^(m)dm is the number of stars 
in the mass interval m to m + dm. 

F or the thin disc, we adopt the IMF of Kroupa et al.l 
(T993E which is similar to that of 'Miller & Scalo" (1979|) and 
IZoccali et al.l (12000) . The primary mass is generated using the 
following formula 



^(m) 



0.035m-'-3 , 
0.019m-22 ^ 



0.08 < m/Mo < 0.5, 
O.5<ot/M0 < 1.0, 
1.0<m/Mo< 100.0. 



To the thick disc and the halo, we also apply a simple IMF 
of power law, 

^(m) oc m-". (35) 

Here, we take or =1.5 for the thick disc and the halo. 

We have adopted a metallicity Z = 0.02 (Population I) for the 
bulge, thin disc and thick disc, and Z = 0.001 for the halo. We 
have also carried out calculations for the thin disc with metallic- 
ity Z - 0.001 in order to see the effect of metallicity on the GW 
signal. 
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Fig. 2. Circular velocity as a function of galactocentric distance 
R from the Galatic model, showing the contribution due to differ- 
ent components, i.e. the bulge, thin disc + thick disc, and halo in- 
cluding dark matter The dashed line indicates the observational 
estimate by Brand & Blitz (1993). The spheroidal component 
due to the interstellar medium was not considered separately. 



(iii) We assume a constant mass- ratio distribution 
(iMazeh et al.lll992HGoldberg & Mazehlll994h . 



n(l/^)= 1,0 < l/q< 1, 



(36) 



(iv) We employ the d istribution of initi al orbital separations 
used by iHanl (Il998l) and iHan et al.l (l2003l) . where they assume 
that all stars are members of binary systems and that the distri- 
bution of separations is constant in log a (a is the separation) for 
wide binaries and falls off smoothly at close separations: 



an(a) = { ^ "0 

asep,ao < a <a\. 



(37) 



where asep ~ 0.070, aq = \QRq, oi = 5.75 x IO^Rq = 0.13 
pc, A; X! 1 .2. This distribution implies that the number of binary 
systems per logarithmic interval is constant. In addition, approx- 
imately 50 per cent of all systems are binary stars with orbital 
periods of less than 100 yr. 

(v) The distribution of eccentricities of binaries follows P,, - 
2e jNelemans et al. 200 Ibl) . 



^^^> 2.5. Rotation curve and local stellar density 



iKroupal (1200 ih gives ^(m) oc m ^-^,0.5 < ot/Mq, see I 



In order to see the influence of the Galactic model on the 
rotation curve of the Milky Way, we plotted Fig. [2] We 
used the Miyamo to-Nagai potential (iMivamoto & Nagailll975l : 
iRevaz et ani2009l) in cylindrical coordinates for calculating the 
circular velocity of the bulge and disc components. For the dark 
matter halo, we adopted the potential of Caldwell & Ostrikea 
(!l98ll). The observational estimate bv .Brand & Blitz. (.1993.) is 
included for comparison. 

From the Galactic model, the total mass of the halo including 
dark matter is 4.5 x 1O"M0 inside a sphere of radius 50kpc. In 
this paper, we only focus on the baryonic mass in the halo which 
is considered to be 1 x IO'Mq. With the SFR adopted here, the 
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baryonic mass in the bulge and disc is at least 2 x 10'"Mq and 
5.5 X 1O"^M0 respectively, implying that our model requires no 
dark matter component in the bulge or the thin disc. 

Combining the Galactic model and the mass of the Galactic 
components, the stellar density in the solar neighbourhood is 
O.O64M0pc"^ of which 6.27 x lO^^Mopc"^ is in the thin disc, 
9.4 X lO^^'Mopc"^ is in the thick disc, and 2.18 x lO^^Mopc"^ 
is in the halo. This i s consistent with the Hipparcos result, 
0.076 + O.OI5M0PC-3 (ICrezeetal.lll998h . The local dark mat- 
ter density in our model is about O.OlMopc"^. 

2.6. Procedure 

In order to compute birth rates, number densities, space, mass 
and orbital distribution of DWDs, we have adopted the follow- 
ing procedure. For each Galactic component g (bulge, thin disc, 
thick disc and halo) having a total mass Mg{t) — L SFR(f')df': 

1 . Calculate a sample distribution of k coeval binaries having a 
total mass m^ and generated by the four Monte-Carlo simu- 
lation parameters m, a, q and e. 

2. Follow the evolution of each primordial binary to establish 
the properties of n(t') DWDs formed from the original sam- 
ple, and the number x{t') of DWDs which merge, where t' 
represents time since formation of each binary. 

3. By combining n(f') and x{t') with the star formation rate, 
compute the birth rate of DWDs v(0 - j^'(«(f')/»ip)xSFR(f- 
f')df'. 

4. Compute the merger rate of DWDs ^(f) = ^ (x(t')/mp) x 
SFR(f - t')dt'. 

5. Evaluate the total number of DWDs N{t) = J^ v(t')-^(t')dt' . 

6. For a given t (e.g. 14 Gyr), use a second Monte Carlo proce- 
dure to generate the orbital properties and spatial distribution 
of N{t) Galactic DWDs by interpolation and extrapolation on 
the n DWDs from the simulation. 

7. Sort the DWDs by orbital frequency. 

8. Calculate the total strain amplitude h^ from the number and 
distance of DWDs in each frequency bin. 

In our population synthesis simulation, we started with k - 
1.20x 10^ primordial binaries (4.0x 10'' in the bulge, 5.0x 10'' in 
the thin disc, 1 .0 x 10'' in the thick disc and 2.0 x 10'' in the halo) 
yielding a total of 7.66 x lO'* DWDs over all four components of 
the Galaxy at the present epoch, t - 14 Gyr 



3. Results 

3.1. Birth rates, local densities and numbers of DWDs 

The results of our simulation are shown in Table |2] where we 
give the birth and merger rates (number per year), local densi- 
ties (pc^^) and total numbers of DWDs. Figures are given for 
each component of the Galaxy and for the Galaxy as a whole. 
Table |2] also shows the numbers detectable if we assume that 
DWDs continue t o be d etectable up to 10** yr after their for- 
mation (Iben et al.l |l997h . At the end of Table |2] we show the 
birthrate for Type la supernovae - 0.0013 yr', assuming that 
all SNe la are formed from merging doubl e carbon-oxygen (C O) 
white dwarfs with total mass > 1.378Mo (iMartin et al.ll2006h . 

The local densities of DWDs in our model are 2.1 x 10 "*, 
4.4 X lO"*", and 4.1 x lO^^pc"^ in the thin disc, thick disc, and 



halo, respectively. The predicted density of halo DWDs is sig- 
nificantly less than the observed density of halo white dwarfs 



= 2 X 10 "^pc^^ (lOppenheimer et aT]l200lh . This could be par- 
tially due to the star formation history, as we only have a single 
star burst of 1 x lO^M© at the formation of the halo. If we as- 
sume 10% of halo WDs are in binary systems (Holberg 20091), 
the observation requires the local density of halo DWDs to be 
a; 2 X lO^^pc"^. This would require the baryonic mass in the 
halo to be 5 X IO'^Mq, 11% of the total mass of the halo, while 
the dark matter takes up the remaining 89% of mass of halo. 

Table |2] also gives figures for the different thin disc models 
discussed above. A low metallicity results in a slight increase of 
the birth rate and n umber of W D binaries for a given IMF, in line 
with the results of Han ( 1998). 

Figure [3] shows the distribution of orbital period Poib as a 
function of primary mass (mi, left panel) and mass ratio (q = 
m2/mi, right panel). 



3.2. Chirp masses 

We plot the number density distribution of DWD chirp masses 
against frequency in Fig.|4] Using Eq.|23]and assuming 15 Gyr 
for the age of universe, we obtain a relation between chirp mass 
and a critical frequency for DWDs with a circular orbit, i.e. e = 
0.0, 

A1 = 3.05xlO"Vc"'*^^- (38) 

Figure |4] shows this relation; DWDs with / > /c will merge 
within 15 Gyr. 

We can understand the distribution of chirp mass from the 
view of stellar evolution. A star only develops a degenerate he- 
lium core to become a helium WD either if its main sequence 
progenitor evolves to the giant branch and loses its hydrogen en- 
velope prior to core helium ignition, or if the hydrogen-envelope 
is removed during main-sequence evolution so that the star fails 
to reach the giant branch and evolves directly from either the 
main sequence or the Hertzsprung gap to the He WD cooling 
track. The mass of the He WD will be < O.5M0, depending 
somewhat on metallicity. 

If the mass of the He core is > 0.5Mo, core helium burn- 
ing will be ignited. Following core-helium exhaustion, the star 
will evolve to become a degenerate carbon-oxygen (CO) or 
oxygen-neon-magnesium (ONeMg) WD with mass in the range 
of =s 0.5 - 1 .44M0 following an asymptotic giant branch, naked 
helium giant or hot subdwarf phase, depending on the chemical 
composition and mass of the progenitor. 

With the population synthesis parameters adopted here. 
Figure |4] shows a principal ridge in chirp mass at A1 ~ O.3M0. 
Figure |6] demonstrates that this ridge corresponds to the He-i-He 
DWDs and represents the largest population of DWDs in our 
model. A smaller and broader ridge is located at M ~ O.65M0 
and con-esponds to the CO+CO DWDs (Fig.©. CO-hHe DWDs 
have a range of chirp masses intermediate between the He-i-He 
and CO+CO DWS (Fig.|6l), making the chirp-mass distribution 
in Fig.|4]appear relatively continuous. 

3.3. The total GW amplitude spectrum from DWDs and 
comparison with LISA 

In order to synthesize the amplitude spectrum for DWDs, we 
choose a LISA integration time T = 3.16x 10^ s (1 yr), so that 
the frequency resolution Af - l/T - 3.17 X 10^^ Hz, which we 
adopt for the size of the frequency bin. Hereafter, GW frequen- 
cies / are given in Hz, unless otherwise stated. 

Some frequency bins in the high-frequency domain will con- 
tain an individual source if the bin is small enough, arising 



S. Yu & C. S. JefFery: Gravitational Waves from Double White Dwarfs 







1.2 

1.0 

0.8 

0.6 

0.4 

0.2 

0.0 

-0.2 

-0.4 

-0.6 

-2.5 -2.0 - 




1.5 -1.0 -0.5 0.0 0.5 1.0 1.5 2.0 2.5 

logPorb (days) 



5 -2.0 -1.5 -1.0 -0.5 0.0 0.5 1.0 1.5 2.0 2.5 

logPorb (days) 



Fig. 3. Distribution of primary masses (left panel), mass ratio q - m2//«i (right panel) and orbital periods of DWDs in our model. 
The contour scale to the right of the figure represents the number of systems in each bin. The present DWD population is plotted, 
but the total number of DWDs (2.76 x 10**) is reduced by a factor of 100. Filled triangles and squares are for known WD-t-WD 
systems. The inset in the left panel shows the distribution of primary mass and orbital periods at logPoib < -1.6 days. Bin sizes are 
Alog(Porb/day) = 0.05, Ami /Mo = A^ = 0.02. 

Table 2. Birth rates, local densities and numbers of DWDs. v = current birth rate, ^ - current merger rate, A^ - total number, A^o = 
the number for DWDs with a detectable component, X - merger number, pLD=local density of DWDs. The last two rows represent 
the current birth rate of supernovae la and the number of resolved DWDs detectable by LISA. Different thin disc models are denoted 
by a: Z = 0.02; b:Z - 0.001. Rates v, 4^are in yr~', pld in pc'^ ■ Note: a DWD is deemed detectable f or 10 ^ yr aft er the second star 
became a WD. This definition foUows Hben et alj d 19971) and the values are given for comparison with lHanI (Il998l) . 



N 
No 

X 

Pld 



Bulge Thin disc 

8.0 X 10-^ "2.1 X 10-^ 
*2.7 X 10-2 
n.7x lO** 
*2.3 X lO** 
"2.1 X 10*^ 
*2.7 X 10*" 
n.4x 10-3 
*1.5x 10-3 
n.lx 10' 
*l.lx 10' 
"2.1 X 10-* 
*2.8 X lO-* 



Thick disc 



Halo 



7.6 X 10' 
8.0 X 10^ 
4.0 X 10-* 
5.4 X 10« 



1.4 X 10-3 
1.2 X 10' 
1.4 X 10^ 

1.1 X 10-* 

1.2 X 10" 
4.4 X 10-* 



1.7 X 10-3 
1.9 X 10' 
1.7 X 10^ 

8.3 X 10-* 
1.2 X lO** 

4.1 X 10-' 



Galaxy 
"3.21 X 10-2 
''3.81 X 10-2 
"2.76 X 10** 
''3.37 X 10** 
"3.21 X 10* 
''3.81 X 10* 
"2.74 X 10-3 
''2.84 X 10-3 
"1.88 X 10' 
''1.88 X 10' 
"2.2 X 10-* 
'^2.9 X 10-* 



SNela 






v = "1.3x 10-3 


resolved 


14800 


"18870 


"A', = 33670 






''22670 


''N, = 37470 



from a rapid decrease in the number of Galactic binaries to- 
wards short e r perio ds. These are the so-called resolved sources. 
lEvans et aP (Il987l) discussed the relation between the detector 
bandwidth A/ and integration time T. They plotted the ampli- 
tude of DWDs with integrations times of 10^ s and 10^ s, and 



found that all binaries with a frequency above 
be resolved if the integration time is > 10^ s. 



3 mHz should 



Figure |5] illustrates the predicted total number of DWDs per 
frequency bin (right panel) and the GW signal (left) these DWDs 
would produce. Some bins in the range -3.0 < log/ < -1.80 
only contain a single system (Resolved DWDs), making up 
-0.012% (33670, Tablellli of the total number of DWDs. 

LISA is designed to be a space-based GW detector, con- 
sisting of 3 satellites flying in formation to form a Michelson 
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Fig. 4. Distribution of the chirp mass and frequency of DWDs. The contour scale to the right of the figure shows the number of 
systems in each bin. The present DWD population is plotted, but the total number of DWDs (2.76 x 10**) is reduced by a factor of 
100. The solid line denotes the boundary for DWDs which will merge within 15 Gyr. The upper-right panel shows the distribution 
of chirp mass at high GW frequency log/ > -3.0 Hz. Bin sizes are Alog/ - 0.05, and AA1 - O.O2M0. 
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interferometer with an arm length of 5x10^ km. Noise arises 
mainly from the laser tracking system (position noise) and para- 
sitic fo rces on the proof m ass of an accelerometer (acceleration 
noise) dLarson et al.ll200C)l) . We can convert the noise signal to an 
equivalent GW signal in frequency space by 



(39) 



where 5n is the total strain noise spectral density, hi is the root 
spectral den s ity an d R is the GW transfer function given by 
iLarson et all (l2000h . 

For a continuous monochromatic source, such as a DWD 
with a circular orbit, which is observed over a time T, the root 
spectral density will ap pear in a Fou rier spectrum as a single 
spectral line in the form dLarson et al. 2000) 



hf 



h 



A/ 



/iVf. 



(40) 



So, for an observation time T - I yr, the root spectral density 
/if = 5.62 X \(fih. 

To demonstrate the detectability of the predicted GW signals 
due to DWDs, Figs. |5]-|9] show the expected LISA sensitivity 
forS/N=land5. 

3.4. The GW signal from various DWD types and cliannels 

Figure |6] shows the contribution to the GW signal from differ- 
ent types of DWD in our model, including those containing two 
helium WDs (He-i-He), a carbon-oxygen WD and a helium WD 
(CO-i-He), and two CO WDs (CO-hCO). DWDs containing at 
least one ONeMg white dwarf are designated "ONeMg". 

Figure |6] shows that CO+CO and ONeMg DWDs have a 
stronger total strain amplitude for log/ < -3.5 (Porh > 0.9 hr), 
up to almost 0.6 dex higher than that of He-i-He and CO-i-He 
DWDs, despite being less numerous. Since the strain amplitude 
for a single DWD of given frequency is proportional to A1^^^ 
(Eq. |7]), the larger chirp masses of the CO+CO and ONeMg 
DWDs dominate the numerical superiority of the He-i-He DWDs. 

CO+CO and ONeMg DWDs are formed in more massive 
progenitor binaries than He+He and CO+He DWD progenitors. 
Such binaries undergo RLOF or CE ejection when the stars are 
physically larger, and hence are in longer period systems and less 
likely to form short-period DWDs. Consequently the slope of the 
number-frequency distribution for CO+CO and ONeMg DWDs 
is steeper than for less-massive DWDs, and the nett contribution 
to the strain amplitude falls away more quickly. 

Consequently, for log/ > -3.5, the signals from all four 
DWD types converge and then at log/ > -3.0, the strain am- 
plitude appears to scale roughly as the average chirp mass for 
the DWD type. The main reason is that the number of systems 
per frequency bin is here very small and generally < 3. The con- 
tribution from ONeMg DWDs remains highest because of the 
chirp-mass effect. However, our model predicts negligibly few 
ONeMg DWDs at log/ > -2.5, and negUgibly few CO+CO 
DWDs at log/ > -2.25 (Porb < 6.65 min), so these contribu- 
tions vanish at these frequencies. 

Fig. |6] shows that CO+He DWDs are the dominant GW 
source at very high frequency (log/ > -2.3). This may be 
understood as a consequence of evolutionary age. In general, 
ONeMg DWDs and CO+CO DWDs coiTespond to larger pro- 
genitor masses and hence a smaller number of progenitors, faster 
evolution and a lower detection probability. On the other hand. 



He WDs have larger radii, so He+He DWDs tend to merge more 
quickly after a double CE phase. 

Figure|7]compares the contribution of different DWD forma- 
tion channels. The stable RLOF+CE and CE+CE channels dom- 
inate the production of DWDs, with the CE+CE channel gener- 
ating 66.7% of all DWDs and dominating the GW signal at all 
frequencies log/ > -5 except in the range -3.0 < log / < -2.3 
where the contribution from the RLOF+CE and CE+CE chan- 
nels may be comparable. The exposed core + CE channel is a 
very minor channel (< 3%) and has almost no influence on the 
GW signal at frequencies log / > -5 (Porb < 55.56 hr) in our 
current model, but this would depend on a tidally enhanced stel- 
lar wind as discussed by Han ( 1998) (see also ^2.2.6l l. 

FigureQalso shows that the highest-frequency DWDs come 
exclusively from the double CE channel. This channel is ex- 
pected to be the most important channel for close DWDs be- 
cause a deep spiral-in phase is necessary to release sufficient or- 
bital energy to eject a common envelope, finally leading to the 
formation of a very close binary. 

In Table [3j we list the total numbers and relative contribu- 
tions of various types of DWD and evolutionary channels, as 
well as their current birth rate and the potential number of re- 
solved systems. We find that CO+He and He+He DWDs would 
be the most significant GW sources at log / > -2.92 (Porh < 
27.7 min), at which frequency resolved sources emerge from our 
sample. We however cannot neglect the population "ONeMg" as 
they also have a considerable resolved number. CE+CE is the 
dominant formation channel for resolved DWDs, producing 3.8 
times more systems than the stable RLOF+CE channel. 

3.5. Tlie GW signal from various Galactic populations 

Figure [8] shows the GW amplitude (right panel) and the number 
of DWDs per frequency bin (left panel) due to each component 
of the Galaxy, including bulge, thin disc, thick disc and halo. 
This figure demonstrates that DWDs in the thin disc and bulge 
generate a strong strain amplitude which should be observable 
by LISA. DWDs in the halo and thick disc might make a sub- 
stantially smaller contribution to the GW amplitude and only at 
log/ < -3. These DWDs will not contribute to the LISA back- 
ground. 

Figure |8] shows that the main difference in the GW signal 
between different populations is at very high frequencies. At 
lower GW frequencies the discrepancy is principally caused by 
the number of DWDs. At very high frequency, the GW signal 
strongly depends on the star formation history. We note that a 
persistent quasi-exponential SFR was applied to the bulge and 
thin disc, but only a single star burst was applied at the begin- 
ning of the thick disc and halo. Figure [8] also shows that there 
is only a slight difference between the signal from the bulge and 
thin disc, despite having a differnt IMF and spatial distribution. 

Figures [1] and |8] show that distance has a only small effect 
on the DWD signal at high frequency. Therefore, if we assume 
d ^ 10 kpc, take logarithms of Equation[8] and combine with the 
relation between orbital period and GW frequency (/ - 2/Porb'- 
circular orbits), we obtain 



log /.. -20.13 + ^ log (^) + ^ log/. 



(41) 



We plot this relation in Fig. [8] as black dashed lines, which 
demonstrates the effect of chirp mass on GW amplitude. Thus at 
high frequency, where individual DWDs can be resolved, chirp 
mass will strongly influence the strength of the GW signal. It 
will be used to compare our results with others in §|4] 
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Fig. 6. The GW signal due to different types of DWD. log/; represents the strain amplitude (left panel). The number distribution 
is shown in the middle panel. The right panel illustrates the distribution of chirp mass due to different DWD types, normalized by 
the total number (2.76x10**) of DWDs in the Galaxy (bin size AM = O.OOIM0). The LISA sensitivity and frequency bins are as in 
Fig.S 
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Fig. 7. The GW signal due to DWDs from different formation channels. The left panel presents the strain amplitudes. The number 
distribution is shown in the right. The LISA sensitivity and frequency bins are as in Fig.|5] Note that the exposed core + CE channel 
is a very minor channel (< 3%) and has almost no influence on the GW signal at frequencies log/ > -5 (Porb < 55.56 hr). 



3.Q. The effect of metallicity 

There is a large number of parameters in our model for synthe- 
sizing the evolution of large numbers of stars and, since every 
synthesis takes a large amount of computing time, too many to 
explore individually in this paper However, by far the most im- 
portant parameter in any stellar evolution model, after mass, is 
metallicity. This will therefore be a crucial factor in determin- 



ing the properties of the D WD population dPols et al.lll998l : lHanl 
[T99ilHurievetani2000l) . 

In order to show the influence of metallicity on the GW sig- 
nal, we performed an additional computation for the thin disc 
using identical model parameters except for the metallicity. We 
have compared metallicities Z = 0.02 and 0.001. 

Figure |9] demonstrates that the shape of the GW spectrum is 
slightly affected by metallicity at low frequency, but is not much 
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Table 3. Birth rates of DWDs with various types and different formation channels, v - current birth rate, A^ = total number, % = 
percentage of the type of DWDs making up the total resolved DWDs in the Galaxy, NResoived - number of resolved DWDs. We take 
the thin disc model with Z - 0.02. 





V 


N 


N/N,,i 


^Resolved 


% 


types of WD binaries 












He+He 


1.34 X 10-2 


1.06 X 10** 


38.41% 


19936 


59.21% 


CO+He 


5.07 X 10-3 


4.11 X 10' 


14.89% 


12852 


38.17% 


CO+CO 


1.15 X 10-2 


1.08 X 10** 


39.13% 


586 


1.74% 


ONeMg 


2.13 X 10-3 


2.09 X 10' 


7.57% 


296 


0.88% 


formation channels 












RLOF+CE 


9.86 X 10-3 


8.27 X 10' 


29.96% 


1061 


3.15% 


CE+CE 


2.12 X 10-2 


1.84 X 10** 


66.67% 


32609 


96.85% 
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9.3 X 10« 


3.37% 
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Fig. 8. The GW signal due to DWDs from different populations (see key). The black dashed lines denote the relation of amplitude 
and frequency fromEq.l4T]at chu-p masses 0.12, 0.22, 0.32, 0.42, 0.52, 0.62, 0.72 and 0.82 Mq (from bottom). The LISA sensitivity 
and frequency bins are as in Fig.|5] 



affected at high frequency. A lower metallicity increases both the 
total number of DWDs, and the number of resolved DWDs. In 
our thin disc model, these numbers are 1.7 X 10^(2.3 x 10'^) and 
18870 (22670) forZ = 0.02(0.001), respectively. 

The difference in the shape of the spectrum arises because at 
lower metallicities, asymptotic-giant branch stars are able to de- 
velop more massive cores, and also achieve higher luminosities 
and larger radii. This leads to a greater number of more mas- 
sive WDs in longer initial period orbits, and also to more DWDs 
surviving. 



Hubble time or equivalent; we adopt lOGyr Considering only 
the 1.7 X 10*^ thin disc DWDs (Table |2]i this fraction represents 
some 94% (or 1.6 x 10^) of the total. 

2) Of the remaining 6% that will evolve into contact within 10 
Gyr, the majority, 95% or 9.6 x 10'' DWDs are formed with 
q> qc and will merge. 

3) The remainder, having q < qc will undergo stable Roche lobe 
overflow. This will decrease their mass ratios further, their orbits 
will widen, and the systems will only stay in contact by further 
GW radiation. These correspond to some of the AM CVn sys- 
temo Their number is a; 5% of those that evolve into contact. 



3. 7. Detached, semi-detached and merged DWDs 

Following formation of a DWD, evolution under the influence of 

GW may follow one of several paths: 

1) If the initial binary is too wide, it will not merge within a 



3 iNelemans et alj ( 12001 ah also discuss AM CVn systems which are 
formed from low-mass helium stars with degenerate companions. 
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Fig. 9. The GW signal obtained for different thin disc models. The red line is for Z =0.02, green for Z =0.001 . The LISA sensitivity 
and frequency bins are as in Fig.|5] 



or 4.8 X 10^ DWDs, or 0.3% of all DWDs formed in the thin 
discfl 



3.8. Resolved sources of GW radiation 



lEvans et all (Il987h suggested that degenerate dwarf binaries 
would be promising sources of GW and pointed out that the de- 
tection of GW is related to the integratio n time. By computing a 
large number of DWDs. lNelemans et al.l OOOlb ) concluded that 
wide DWDs dominate the GW signal at log/ < -3.4, and that 
AMCVn systems could be resolved at log/ > -2.8, but would 
produce a lower signal due to their small number. Subsequently, 
iNelemans et al.1 (|2004|) showed that AM CVn systems would be 
good candidates for LISA and that k 11000 AMCVn systems 
could be resolved. By "resolved", we mean that, over a certain 
observation time T, a frequency bin (A/ =1/7") contains only 
one binary. 

In general, our methods and results are similar One differ- 
ence is that, by including their contribution, we confirm that 



■* These figures correspond to model I of Nelemans et alj ( 12001 ah . 
in which there is no tidal c oupling between the acc retor spin and the 
orbital angular momentum. INelemans et al.l ([2p01a) find that with ef- 
fective tidal coupling, the number of AM CVn systems may increase by 
two orders of magnitude (their model II). It is not the intention of this 
paper to discuss AM CVn systems in detail. 



DWDs from the halo and thick disk would make little contribu- 
tion to the LISA signal. A second difference is that, by adopting 
a different IMF, star formation rate and total mass for the bulge 
and the thin disc, we obtain an increase in the total number of 
compact DWDs at high frequencies, so that the number of re- 
solved sources increases slightly. 

Our results confirm that all DWDs with frequency log / > 
-2.25 would be resolved with a LISA integration time equal 
to 1 yr In §3.3, we noted that resolved systems with log/ > 
-3.0 would account for 0.012% of the total number of DWDs, 
of which 29% will become semi-detached. These numbers 
agree reasonably with previous results, although Nelemans et aP 
(l2001bi 120041) found ^ 50% to be semi-detached. 

93.41 points out that most resolved DWDs should have 
formed from the CE-f CE channel and should be dominated by 
CO-nHe and Hen-He DWDs. They should be the most significant 
resolved sources for LISA due to their number and strain ampli- 
tude and would dominate the LISA GW signal at log / > -2.25, 
(^orb < 6.65 min). It may also be possible to infer their type 
(CO-nHe or Hen-He) from their chirp masses (Eq.l42li by deter- 
mining f orb and Poib, possibly from optical observations. 

CO-hCO and ONeMg DWDs have a stronger strain ampli- 
tude at log/ < -3.12 (Porb > 0.732hr), resulting from their 
larger chirp masses. For frequencies -2.25 > log/ > -3.12, the 
signal from all four types of DWD would overlap. 
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Table 4. Main parameters and their value in our simulation. See Table [T] for the Galactic structure parameters. See 
explanation of the parameters. 
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Table 5. Estimated errors in DWD birthrates for the thin disk 



p 


(P) 


(Sp) 


dv/dp 


Svp/v 


7 


1.5 


0.1 


3.4 X 10-* 


0.016 


logZ 


-1.7 


0.5 


4.3 X 10-3 


0.102 


rGyr-' 


14 


1 


0.1 X 10-3 


0.005 


aiMF 


-2.2 


0.1 


8.6 X 10-" 


0.041 


flSFR 


11 


1 


1.0 X 10-3 


0.047 



3.9. Error Analysis 

Modelling stellar populations and hence estimating the numbers, 
properties and distributions of stars in the Galaxy is subject to 
several sources of error (or uncertainty). These may naively be 
divided into statistical errors arising from the implementation 
of the Monte Carlo process, and systematic errors arising from 
choices made for the model parameters (of which over thirty 
may be identified in § 2; principal values are indicated in Table 
II. 

The Monte Carlo procedures involve two stages, the first of 
which creates a sample of J ndt' - 1.84 x 10'* DWDs in the thin 
disc for example, on which subsequent samples are interpolated. 
Assuming Poisson statistics, the error (Icr) in this first number 

is a: J J ndt' a: 1.4 X 10^ or a 0.8%. This represents a lower 

limit on the eiTor in subsequent quantities, including total birth 
rates, merger rates and numbers of DWDs. The statistical error 
on the number distributions and strain amplitude spectra will be 
at least as large as this, increasing as (A^(/)ii/)-"-^ towards higher 
frequencies. 

In order to investigate some of the systematic errors, we have 
caiTied out reduced simulations for the thin disc population with 
an initial sample of 10^ binaries and by varying each of the five 
most important parameters (Table |5]). Whilst we omit other pa- 
rameters from this exploration (B, for example) to save com puta- 
tion t ime, their systematics have been discussed by others (iHanl 
[T998h . 

Systematic errors can be estimated by considering the gradi- 
ent dv/dp of the birthrate v with respect to a model parameter p, 
and multiplying by an expectation value (dp) for the uncertainty 
in p. Considering y, metallicity (here expressed as logZ), age 
t, IMF, and SFR to be representative of inputs with large uncer- 
tainty, fractional eiTors 6vp/v = dv/dp x {6p}/v are indicated in 
Table |5] To simplify, we have here assumed the disk IMF takes 
the simpler form ^(m) » m" and the disk SFR takes the form 



fl(exp(f - t())/T) + b(t - to), so that they may be characterised by 
ffiMF and flsFR respectively; t is constraine d by the current SFR . 
(foi Mp) ^ . 1 is su ggested by comparing iKroupa et al.l (Il993h 
with iKroupal (1200 li) . while (^cisfr) ~ 1 is taken f r om th e uncer- 
tainty in the thin disc mass given bv lKlvpin et alJ (l2002h . 

The dominant uncertainties in v arise from the metallicity 
approximation Z = 0.02, at least in the thin disc, and the SFR. 
Again making a naive approximation that all p are independent, 
adding the error contributions quadratically gives a total relative 
error (Icr) from these sources 6v/v ~ 0.34. 

These numbers are also useful for comparing results between 
different models by substituting the uncertainty dp with the dif- 
ference between adopted values Ap. 

In our model the Galactic D WD birthrate v - ^ ~ 3.2 x 
10-2 yr' (Table |2]l. JHai (Il998h found a similar value (3 x 
10-2 yr'). For the thin disc alone we find th e DWD birthrate 
to be 2.1 X 10-2 yr' while iNelemans et alJ (feOOlb ) obtained 
2.5 X 10~2 yr-'. The s l ight dif ference between our model and 
that of INelemans et al.l (1200 Ibl) can be understood primarily in 
terms of the difference AasFR = 3 and AaiMp =0.3 

The birth rate of supernovae la (0 . 0013 yr-') in our model 
is sUghtly sm aller than that of ES (Il998l) (0.003 yr-') and 
INelemans et al . (2001b) (0.002 yr"'), and is probably also a con- 
sequence of the IMF and SFR. The fact that we set a limit to the 
Eddin gton accretion rate (Eq. l25l l restricts the SN la rate (iLiviol 
120001) . 

Another effect on the birth rate and number arises from 
tidally-enhanced stellar winds and / or wind-driven mass^ trans- 
fer, which is related to a coeflicient B ( g2.2.2l l. lHml(ll998l) finds 
that the birth rate and number vary slightly with B. We adopt 
B - 1000 following Han (1998), but much large r values might 
be justified; with B = 10000 dTout & Hallll 19911) . the mass-loss 
rate could be 150 times larger than the Reimers' rate when the 
star nearly fills its Roche lobe. 



4. Discussion 



In addition to work by INelemans et al.1 (12001 hi |2004 already 
discussed, the question of GW radiation from DWDs and their 
contribu tion to t he LISA signal has al s o bee n addressed by 
Webbink & Had (1 19981 
(2000), Willems etail 
(l2009l) . amongst others. 
In our model, the 



iHils & Beiidd (l2000l) . iHiscocketal 
(120071) . iLiul (120091) and iRuiteretal 



GW 



Galactic DWD population would be 
10-24 Hz'/2 for log/ > -4.0, and 



strain amplitude from the total 
2.51 X 10-22 - 6.31 X 
5.01 X 10-22 - 5.01 X 
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10"^^ Hz'^^ for log / < -4.0. The larger strain at lower frequen- 
cies is primarily due to the larger number of DWDs per unit fre- 
quency. 

At log/ > -3, the GW signal from DWDs in the bulge and 
thin disc combined equals the total Galactic GW signal. The 
thick disc and halo make no contribution (Fig. [8]l because all 
DWDs that formed in these components that will m erge have al- 
ready done so. At lower frequencies (log/ < -3). iRuiter et alj 
(l2009l) found that the halo signal at log/ = -3, for example, 
hf = 10-20 Hzi/2(= iog/iHz-i/2 ^ -23.75 ) is a factor ^ 10 

lower than that of the disc-nbulge, where hf - 10-'^Hz'^2(= 
log/iHz"'/2 = -22.75). They estimated there to be 1.5 x 10'' 
halo white dwarfs with a total mass of lO'^M© in a halo with 
Z = 0.0001 and age 13 Gyr. At log/ = -3.3 we have < 
log/i > ^ -23.3, roughly a factor three stronger than^ Ruiter et alJ 
(120091) . but in a halo with Z = 0.001 and age 14 Gyr. 

Another differe nce between our model and the model of 
iRuiter et al] (l2009f) arises from the chirp masses (cf. § 13.51 ). 
Figure [8] shows the strain-amplitude relation as a function of 
chirp masses; since there are only one or two systems in the ma- 
jority of high frequency bins, it is possible to infer that the chirp 
masses of most DWDs in our model are greater than 0.1 Mq. In 
comparison, the chirp masses of high-frequency DWDs in the 
model of Ruiter et al. (2009) appear to be less than 0.1 Mq. The 
different distribution of chirpmass in the two models may be due 
to differences in Z and age adopted in each case. Ou r results 
more closely resemble those of Nele mans et alj (l2001bh . 

The GW radiation from AM C Vn binaries and helium cata- 
clysmic variables was calculated bv lHils & Bendeij ( 12000 ), who 
found that these two populations only provide a slight enhance- 
ment on the confusion noise for LISA below 3 mHz and even 
no increase at higher frequency. Nelemans et al. (2001b) and 
IWebbink & H an ( 1998) showed a confused background GW sig- 
nal due to DWDs which is similar to ours at low frequency, but 
different at high frequency. 

iHiscock et al.l (12000 ) found that the background signal from 
halo WD binaries could be five times stronger than the ex- 
pected contribution from Galactic disk binaries by assuming that 
the fraction of WDs in binaries is the same in the halo as in 
the disk. This suggested the possibility of a GW signal from 
halo DWDs although probably the assumption overestimated the 
num ber of DWDs in halo, which would justify further investiga- 
tion. I^Uemsietal] (l2007h considered GW radiation from eccen- 
tric DWDs formed through interactions in globular clusters and 
found LISA could provide unique dynamical identifications of 
these systems. This requires further investigation. 

Since the change of GW frequency (or orbital period) / can 
be measured by electromagnetic observations or LISA, A4 is 
given by 



/ = 5.8 X 1O"'(M/M0)'/Y 



5/3 f 11/3 



(42) 



Hence, for a small number of high-frequency systems likely to 
be detected by LISA, it will be possible to determine Al di- 
rectly. By knowing M and P, we can calculate the distance d 
from Eq. |7] which is vital for constraining the distribution of 
DWDs and obtaining a better understanding of the structure of 
the Galaxy. 

In this paper, we have not considered the BH in the center 
of the bulge or low-mass BHs elsewhere in the Galaxy. More 
attention should be paid to these in future work, as they could 
greatl y interfere with the abil ity of LISA to detect resolved sys- 






tems jNelemans et al.ll2001bl) . 



5. Conclusion 

A self-consistent numerical simulation of the gravitational wave 
signal from the entire population of double white dwarf bina- 
ries (DWDs) in the Galaxy has been carried out. We have used 
a population synthesis approach to follow the stellar evolution, 
and a comprehensive Galactic model in which we suppose that 
the Galaxy is composed of a bulge, thin disc, thick disc, and 
halo. This model demonstrates a significant contribution to the 
GW signal due to DWDs. 

We have discussed the birth rates, local densities and total 
numbers of DWDs in each fraction of the Galaxy, as well as of 
supernovae la, noting that these are affected by the particular 
choice of IMF, SFR, metallicity and age adopted in our model 
(TableH. 

The GW signals from different components of the Galaxy 
have been computed, as well as the contribution from various 
types of DWD and from DWDs formed by different formation 
channels. The corresponding birth rates and numbers have been 
given (Table [3]l. 

Formation channels CE-hCE and stable RLOF+CE play the 
most important role to produce DWD populations detectable by 
LISA at log/ > -4.5 (Poib < 17.57 h). DWDs with the shortest 
orbital periods in our model come from the CE+CE channel. 
The Exposed Core-nCE channel is an almost negligible channel 
for the formation of DWDs. 

We find that CO-nHe and Hen-He white dwarf pairs dom- 
inate the GW signal at high frequency (log/ > -2.3), while 
CO-i-CO and ONeMg DWDs make the major contribution at 
log/ < -2.3. 

We find that ^ 33 000 DWDs could be resolved by LISA. 
Most would have formed through the CE-hCE channel, and the 
model favours the majority being CO+He and He-t-He DWDs. 
Ground-based observations might determine which by measur- 
ing their chirp masses. 
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